Scalar field induced oscillations of neutron stars and gravitational collapse 
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We study the interaction of massless scalar fields with self-gravitating neutron stars by means of fully 
dynamic numerical simulations of the Einstein-Klein-Gordon perfect fluid system. Our investigation 
is restricted to spherical symmetry and the neutron stars are approximated by relativistic polytropes. 
Studying the nonlinear dynamics of isolated neutron stars is very effectively performed within the 
characteristic formulation of general relativity, in which the spacetime is foliated by a family of 
outgoing light cones. We are able to compactify the entire spacetime on a computational grid and 
simultaneously impose natural radiative boundary conditions and extract accurate radiative signals. 
We study the transfer of energy from the scalar field to the fluid star. We find, in particular, that 
depending on the compactness of the neutron star model, the scalar wave forces the neutron star 
either to oscillate in its radial modes of pulsation or to undergo gravitational collapse to a black hole 
on a dynamical timescale. The radiative signal, read off at future null infinity, shows quasi-normal 
oscillations before the setting of a late time power-law tail. 

PACS number(s):04.25.Dm, 04.40.-b, 95.30. Lz, 04.40.Dg, 97.60.Lf 



I. INTRODUCTION 

Obtaining reliable estimates for gravitational wave sig- 
nals emitted from the stellar collapse of massive stars is 
one of the key motivations for numerical relativity. The 
strength of actual astrophysical sources is still under in- 
vestigation, as the relevance of such sources for the first 
and second generation of interferometric detectors de- 
pends on the details of signal amplitude and frequency, 
but also on the occurrence rates (for a current view, 
see |0). In any case, relativistic collapse is a funda- 
mental physical process, and the development of rele- 
vant computational procedures has been a long steady 
process over the past three decades (see, e.g. B-0]; see 
also II and references therein). In contrast to a Newto- 
nian approximation, where the computational problem is 
well defined and attention can be devoted to astrophys- 
ical details |5|||, there is no consensus as to what is the 
optimal, or at least adequate, framework for developing 
relativistic simulations. 

We will use here the so-called characteristic formula- 
tion of general relativity Jt],||]. The formalism has been 
developed specifically for addressing ambiguities concern- 
ing gravitational radiation and is well adapted to handle 
the propagation of signals. It allows for spacetime com- 
pactification, which avoids problems due to the artificial 
reflection of the fields at outer boundaries. In addition, 
it allows for the extraction of physically relevant global 
quantities like the News function and the Bondi mass. 
Nevertheless, the computation of the dynamics of sources 
of signals is a separate issue altogether, which has not 
been addressed in this early work. The framework for 
computing a complete spacetime within the characteris- 
tic approach has been laid out in Q and more explicitly 
in (lOfl . The translation of this framework into a compu- 



tational tool for vacuum spacetimes, in any dimensions, 
has been a largely successful process (see and refer- 
ences therein). 

There are some noteworthy issues. Firstly, the domain 
of applicability is limited to configurations in which the 
set of light-cones that forms the basis of the coordinate 
system does not fold itself into caustics. This puts a 
limit, for example, on the type of binary systems that 
can be studied. Secondly, to obtain a complete space- 
time one must include in the computational domain the 
vertex of the light-cones. This involves regularity condi- 
tions (and for explicit integration methods, severe time- 
step restrictions), which at present have been resolved 
only up to axisymmetric configurations p2[ . Neverthe- 
less, the approach has remarkable economy and stability, 
which makes it a good candidate for studies of isolated 
relativistic objects emitting gravitational radiation. For 
the simulation of realistic sources, one would have to in- 
clude suitable matter models. It was demonstrated in |l3} | 
that the modern techniques of high-resolution shock- 
capturing schemes for solving the equations of relativistic 
fluid dynamics can be effectively integrated within this 
framework. A separate study implementing the approach 
focused on the gravitational radiation properties of an 
accreting black hole |14|. 

This work is the first example of the use of charac- 
teristic numerical relativity for the study of dynamical 
neutron star spacetimes, collapse and radiative signals. 
The present numerical study is performed in spherical 
symmetry and uses a self-gravitating, massless, scalar 
field. The scalar field serves as a simple matter model 
which mimics gravitational waves. It has been used fre- 
quently to study global properties of spacetimes, black 
hole formation and radiative signals. For the latter, 
the work focused on the interaction of scalar waves and 
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black holes (e.g. JTslJlql ) and especially on the emergence 
of power-law tails ||17|-|20||, which arise from late time 
backscattering of the scalar field at the exterior spacetime 
geometry [^lj . There are few studies of the interaction of 
scalar fields with fluid stellar objects. Recently [^2| an- 
alyzed the scattering of scalar fields off boson stars and 
the emergence of critical solutions for this setup. They 
found that the scalar field can either make the boson star 
collapse to a black hole or to disperse its mass to infinity. 

Time-dependent simulations of scattering of gravita- 
tional waves packets with neutron stars, as a means of 
computing the frequency spectrum of the neutron star 
(see, e.g. pS| for a recent review), have been studied by 
Allen et al"|24| for polytropic equations of state (EoS), 
and by Ruoff |25|] for more realistic EoS. Such simula- 
tions were performed using linear perturbation codes. 
Pavlidou et al. |2(| studied the radiative falloff of scalar 
fields in neutron star spacetimes, using (idealized) ana- 
lytic, constant density neutron star models and assuming 
stationarity for the fluid and the geometry. 

In the present work we investigate the nonlinear dy- 
namics of neutron stars interacting with scalar fields. 
We are especially interested in the following questions: 
How does a stable neutron star react when it interacts 
with the scalar field? Can the scalar field force the star 
to undergo gravitational collapse? What is the result 
of the interaction on the scalar field? In order to an- 
swer these questions we perform spherically symmetric, 
coupled evolutions of the Einstein-Klein-Gordon perfect 
fluid system. We study the reflection of finite scalar wave 
packets off neutron stars for a series of stars parametrized 
by the central density for a given polytropic EoS with 
polytropic index n = 1. Our study focuses on the fate 
of the system during the interaction, both on the gener- 
ation of nonlinear oscillations in the neutron star and on 
the gravitational collapse of the neutron star to a black 
hole. 

The paper is organized as follows: Section [II] describes 
the basic mathematical foundations of our approach. In 
Section [II, we briefly discuss the numerical techniques 



and the implementation used in the simulations. In Sec- 
tion 



of our setup follows the lines of the Tamburino-Winicour- 
formalism in particular as it is applied in regular 
spacetimes, where the foliation of li ght- cones emanates 
from a freely falling central observer (Hj] . 



A. Einstein equations 

By adopting the Bondi-Sachs |?]|| form of the metric 
element in spherical symmetry, 



ds 2 = 



5 2/ V 



du 2 - 2e 2(3 dudr + r 2 {d6 2 + sin(9 2 d</> 2 ), (1) 



the spacetime geometry is completely described by the 
two functions (3{u,r) and V(u,r). 

A sufficient set of Einstein equations for obtaining the 
spacetime development are grouped as 



G 1 — t\H. , 
G uu \t = nT uu \r, 



(2) 
(3) 
(4) 



IV, we present several numerical tests of our code, angular coordinates x 2 



where the u coordinate is defined by the level surfaces 
of a null scalar (i.e., a scalar u satisfying V^iiV^w = 0). 
The r coordinate is chosen to make the spheres of rota- 
tional symmetry have area Airr 2 . The x 2 ,x 3 coordinates 
in this geometry are simply taken to be the angular co- 
ordinates (0, <fi) propagated along the generators of the 
null hypersurface, i.e., they parameterize the different 
light rays on the null cone. With our choice of units 
the constant k is simply n — 87r. The first two Einstein 
equations, Eqs. (g) and (||), contain only radial deriva- 
tives and are to be integrated along each null surface. 
The last equation (^) is a conservation condition, satis- 
fied on the vertex of the null cones T due to the regu- 
larity conditions. We choose T to be a timelike geodesic 
which coincides with the origin of a neutron star at r = 0. 
Equation (0) may be substituted for by the equivalent ex- 
pression g R AB = 8Trg AB (T A B - SUbT/2), where 
is the Ricci tensor and the indices (A, B) run over the 



aimed to assess the correct implementation of its different 
components, the hydrodynamic evolution, the scalar field 
evolution and the metric solver. Section ^ describes the 
actual numerical investigation on the interaction between 
the neutron stars and the scalar field. Finally, Section VI 
summarizes our findings. Throughout the paper we use 
geometrized units G = c = 1 and further assume that 
Mq = 1. Greek indices run from to 3. 



II. MATHEMATICAL FRAMEWORK 

We consider a general spherically symmetric spacetime 
with a two component stress energy tensor of a perfect 



Using the line element and Eqs. (|2j) and (|3j) the (3 and 
V hypersurface equations are given by 



{3. r — 2-KrT rr , 

V r = e 2 "(l~A7rr 2 (g AB T AB -T)) 



(5) 
(G) 



The comma in the above equations indicates, as 
usual, partial differentiation. Boundary conditions for 
(/3(u)r, V(u)r) needed for the radial integrations are pro- 
vided by imposing regularity at the origin, where the co- 
ordinate system is assumed to be a local Fermi system, 
leading to 



fluid and a scalar field, T» v = Tt v 



The geometry 



f3 = 0(r 2 ), 
V = r + 0(r 3 ). 



(7) 
(8) 
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By imposing such condition at the origin, the lapse of 
coordinate time du is related to the corresponding lapse 
of "retarded time" dr measured by distant observers at 
r — > oo as 



dr 



e 2H du, 



where 



H = lim (3. 



B. Scalar Field equations 



(9) 



(10) 



The dynamics of a scalar field $ is governed by the 
minimally coupled Klein- Gordon equation in spherical 
symmetry, 



V M V$ = 0, 



(11) 



where V M is the covariant derivative. This is the equation 
of motion for the Lagrangian 



(12) 



with a corresponding stress energy tensor given by 

T£ v = V M $V$ + Lg^ . (13) 

Using a characteristic foliation the scalar wave equa- 
tion, Eq. (O), takes the form 



2(r$,„), r = -(rV*, r ) 



(14) 



In terms of the intrinsic 2-metric of the (u, r) sub- 
manifold, 



V 



r]cDdx c dx D = -e 2/3 du(—du + 2dr) 



(15) 



where the indices (C,D) run over the coordinates (u, r), 
Eq. (PH) reduces to 



r \r J _ r 



(16) 



where g — r$ and cK 2 ) is the D'Alembertian operator 
associated with tjcd- 



where Tp V is the stress energy tensor of a perfect fluid 



T£" = p/i/u" + pg^ . 



(19) 



All quantities in the above expression have their usual 
meanings: p is the mass density, h = \ + e -\- p/ p is the 
specific enthalpy, e is the specific internal energy and p is 
the pressure of the fluid. Moreover, u M is the four-velocity 
which satisfies the normalization condition g^u^u" = 
-1. 

Following jla], after introducing the definitions D = 
pu°, S r = Tjfand E = T F °, the fluid equations can be 
cast into a first-order flux-conservative, hyperbolic sys- 
tem for the state- vector U = (D, S r ,E) : 



d m + f: 



E.u + F r r — 



(\n^—g), u D 
-{\n^g), r F r \ 

-{hx^g), u E 
-(lnV^.r^-T" TP, 



(20) 
(21) 
(22) 



where \J—g = r 2 sin 9e 2 @ is the four dimensional volume 
element and T" are the Christoffel symbols. The precise 
form of the vector of fluxes F can be obtained by using 
Eqs. @-(|U) (see also The explicit relations be- 

tween the primitive variables w = (p, e, u r ) and the con- 
served variables U = (D, S r ,E), for a perfect fluid EoS, 
p = (r — l)pe, where T is the adiabatic index of the fluid, 
are given in |H| . 

With the above definitions, the metric equations (||)- 
(^) read, for the combined stress energy tensor of a fluid- 
scalar field system, 



P, r = 27rr(ph(u r ) 2 + ($ r f) : 



V, 



(1 — 4irr 2 (ph — 2p)) . 



(23) 
(24) 



Following p7[ we express the hydrodynamic quantities 
on the right-hand side of Eqs. (|23|)-(|2^) solely in terms of 
the conserved hydrodynamic quantities U. This avoids 
additional iterations when using explicit algorithms to 
solve these ordinary differential equations. 

In summary, the initial value problem consists of equa- 
tions (g), p|), (|^||), (§|,|24|), the scalar field initial 
data $(r, uq) and initial and boundary data for the fluid 
variables {p,e,u r ) on the initial slice Eo (at time uq). 
Those equations and initial data are sufficient for ob- 
taining a global solution to the problem. 



C. Hydrodynamic equations 

The evolution of the fluid is determined by the local 
conservation laws of stress energy and density current 

V M T^ = 0, (17) 
V„(pu")=0, (18) 



D. Global quantities 

Making use of the characteristic formulation of general 
relativity and covering the infinite range of the radial co- 
ordinate with a finite grid allows us to refer to some global 
quantities of the spacetime such as the Bondi mass and 
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the news function. Apart from their physical relevance, 
these quantities can be used in global tests of our numer- 
ical evolutions, as we will show below. 

Instead of extracting the Bondi mass directly at future 
null infinity we use the expression 



M = 4tt / r 2 e- 2(3 T r i,dr 



(25) 



for the Bondi mass at time u in our numerical implemen- 
tation. Similarly, the news can be rewritten as 



1 f°° V 

N = -e- 2H / -<5>, r dr. 

2 Jo r 



(26) 



With these definitions, global energy conservation can be 
established, 

M{u) - M (0) = / -AirN{u) 2 e 2H(il) du. (27) 
Jo 



III. NUMERICAL IMPLEMENTATION 

In order to study the interaction of the scalar field and 
the neutron star in a global spacetime we use nonequidis- 
tant grids for the radial coordinate r. Furthermore, to 
avoid dealing with complicated stencils in the numeri- 
cal implementation, we make use of the following pro- 
cedure, generalizing previous implementations in charac- 
teristic numerical relativity (see, e.g. |fT2| ): Starting with 
an equidistant grid in the coordinate x € [0, 1], we allow 
for a general coordinate transformation r — r(x). Using 
the chain rule we rewrite the partial derivatives appear- 
ing in the above equation with 



(),r=(), 



dx 
dr ' 



(28) 



thus effectively rewriting all our equations in the coordi- 
nate x. Unless otherwise stated wc use the relation 



15x 



1 



(29) 



for all computations presented in this work. Using such a 
coordinate transformation the repartition of grid points 
in the coordinate r is almost equidistant for small radii 
and gets infinitely sparse for x — ► 1, which corresponds 
to future null infinity J + . 

We use a second order Runge Kutta method to solve 
the metric equations (|2^) and (p4|). To determine the 
equilibrium models of our relativistic stars, we also use 
the Runge Kutta method to solve the so-called Tolman- 
Oppenheimer-Volkoff equations, formulated on a null hy- 
persurface as in Q. 

The integration of the evolution equation for the scalar 
field, Eq. Jl6| ) (or equivalently Eq. (|l~4|)), proceeds with 
the specification of initial data $(uo ; r) 011 ^ ne initial null 



cone uq- For the characteristic evolution we have used 
two different algorithms, both marching from the origin 
to the exterior. 

The first procedure is based upon the construction of a 
null parallelogram built up from incoming and outgoing 
radial characteristics |2^]. In this procedure one needs 
first to determine the right hand side of Eq. (|l^) at the 
center of the parallelogram to the desired order of ac- 
curacy. Then, an integral relation between this source 
term and the values of g at the four corners of the paral- 
lelogram - which do not necessary have to coincide with 
grid points - has to be employed in order to compute the 
scalar field at that corner of the parallelogram lying next 
to the grid point which is to be updated. Suitable inter- 
polations then give the scalar field at the new grid point 
to second order accuracy. 

The second alternative procedure we have imple- 
mented to solve the scalar field equation is based on a 
direct discretization of Eq. ([TJ) using a second order, fi- 
nite difference, nondissipative algorithm discussed in [^9| . 

Due to the stencils of both algorithms, we cannot use 
them at the origin, where a regular behavior of the scalar 
field as $ = a + br + cr 2 is assumed. The linear term in- 
troduces a kink at the origin, but it is necessary in our 
foliation - as can be seen from the analytic solution for 
the wave equation in Minkowski space consisting of an 
ingoing and outgoing wave. Note that the scalar field 
enters the metric only through Eq. (|23|), thus respecting 
the regularity conditions given by Eqs. (g,§). Substitut- 
ing this ansatz into Eq. (|1J) and grouping those terms 
with the same powers of r we find that a M = fe, b u = 1.5c. 
Extracting the coefficients a, b and c on the null cone uq, 
we update a and b to obtain the scalar field at the first 
two grid points of the new hypersurface, which then al- 
lows us to start the marching procedure along the null 
hypersurface with either of the two algorithms described 
above. 

By experimenting with both algorithms, we found 
that, on the one hand, the scheme based upon a direct 
discretization of the wave equation is more accurate in 
the long-term behavior in the interior of the numerical 
domain. This was relevant to resolve the late time fall-off 
behavior of the scalar field, as we describe below in Sec- 
tion |y|. On the other hand, the algorithm based upon the 
null parallelogram is, however, superior close to future 
null infinity, where we regularized the equations follow- 
ing the work of [Q. Therefore, for the results presented 
in this work we have used a "hybrid algorithm" , in which 
a direct discretization of Eq. ( p^ ) is used in the interior 
of the computational domain and the parallelogram al- 
gorithm is used close to future null infinity . 

Concerning the numerical integration of the system 
of hydrodynamic equations, its hyperbolic mathematical 
character allows for a solution procedure based on the 
computation of (local) Riemann problems at each cell- 
interface of the numerical grid. At cell i the state-vector 
U is updated in time (from u n to u n+1 ) using a conser- 
vative algorithm 
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FIG. 1. Stability curve for neutron star models with the 
polytropic equation of state p = Kp T , with K — 100 and 
r = 2. Models lying to the left of the maximum of the curve, 
at about p c — 3.2 x 10 -3 , are stable against gravitational 
collapse. The circles are calculated with our initial data solver 
and are connected by straight lines. We use units in which 
G — c — Mq = 1. 

n „+i = n „ _ ^(F m/2 - F(_ 1/a ) + AuS, , (30) 

where the numerical fluxes, F, are evaluated at the cell in- 
terfaces according to some particular flux-formula which 
makes explicit use of the full spectral decomposition of 
the system. For our particular formulation of the hydro- 
dynamic equations such characteristic information was 
presented in |p^| . 

In more precise terms the hydrodynamics solver of our 
code uses a second order Godunov-type algorithm, based 
on piecewise linear reconstruction procedures at each 
cell-interface |}0| and the HLLE approximate Riemann 
solver pl]j32|| . General information on such schemes in 
relativistic hydrodynamics can be found, e.g. in p3| ] and 
references therein. 



IV. CODE TESTS 

We now present representative results obtained in the 
process of code calibration. The assessment of the numer- 
ical implementation is provided by comparing to previous 
results and by checking global energy conservation tests. 
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FIG. 2. Time evolution of the radial velocity u x at half ra- 
dius of the star. The neutron star model has a central density 
p c = 1.5 x 10~ 3 . The oscillations are essentially undamped 
for the evolution shown, which reflects the small viscosity of 
the hydrodynamic schemes employed. 

A. Null cone evolutions of self-gravitating, stable 
neutron stars 

As a first step to validate our numerical code we start 
by studying its ability to keep the equilibrium of initially 
stable neutron star models. For this purpose we perform 
long-term simulations of such initial data and analyze 
the stability of the code. Furthermore, we use these evo- 
lutions to compute the frequencies of the radial modes 
of pulsation. We compare the frequencies obtained with 
our nonlinear code to results of linear evolutions from 
perturbation theory. 

For all simulations presented in this paper the neutron 
star models are approximated by a polytropic EoS, de- 
fined by p — K p T , with polytropic constant K = 100 and 
adiabatic exponent r=l + l/n = 2. Hence, the index of 
the polytrope is n — 1. For the simulations presented in 
this section we choose two different models with central 
density p c = 1.5 x 10~ 3 and p c — 2.8 x 10~ 3 (recall that 
we are using units in which Mq = 1). Both models are 
located in the stable branch of the central density - total 
mass - diagram (see Fig. |l|). 

When evolving these neutron star models with our nu- 
merical code we are able to maintain them in stable equi- 
librium for thousands of light-crossing times of the star 
without any sign of numerical instabilities. 

To validate the code further we computed the frequen- 
cies of the radial modes of pulsation. For this aim we 
have to allow the star to (radially) contract and expand 
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FIG. 3. Time evolution of the radial velocity u x at half 
radius of the star. The neutron star model has a central 
density p c = 2.8 x 10~ 3 . 

during the evolution. To this end, following p3 (see 
also |p5| ), we surround the star with a few zones rep- 
resenting an artificial "atmosphere" filling an otherwise 
vacuum region. The density in this atmosphere is set 
to sufficiently small values such that its presence does 
not affect the dynamics of the system. The typical val- 
ues we choose are 10 — 10 -8 times the central den- 
sity of the star. Furthermore, to avoid any numerical 
problems due to (shock) heating in the atmosphere (the 
fluid in those zones is not in equilibrium and, therefore, 
it will collapse/accrete onto the neutron star), we follow 
the recipe described in |34| and enforce adiabatic evo- 
lution (by using the polytropic EoS) in the atmosphere 
and in the outer layers of the neutron star (comprising 
the outermost 10 grid points). After each time step, if the 
density has fallen below 1.5 times the density of the at- 
mosphere, the hydrodynamic quantities are reset to their 
atmosphere values. The innermost location where this 
procedure is done defines the radius of the star. As de- 
scribed in more detail in the next section the above values 
of the atmosphere density are small enough to guaran- 
tee conservation of energy despite the artificial resetting 
procedure. 

When evolving our stellar models in time, we find small 
deviations around the equilibrium values due to the dis- 
cretization errors. As a result, the stars oscillate in their 
fundamental radial modes. Figures ^| and ^ show the ra- 
dial velocity at half stellar radius for the above models 
as a function of the retarded time measured by distant 
observers. These simulations were performed with a grid 
of 800 zones covering the complete radial domain. This 



FIG. 4. Fourier transform of the time evolution shown 
in Fig. The peaks in the Fourier transform indicate the 
mode frequencies of the fundamental radial mode (around 
/ = 1.4 kHz) and the first two harmonics. The neutron star 
model has a central density p c = 1.5 x 10 -3 . The dashed 
vertical lines indicate the corresponding frequencies obtained 
with a perturbative (linear) code. The units in the y-axis are 
arbitrary. 

amounts in using about half of the grid in resolving the 
neutron star (We choose this resolution here to allow for 
comparisons with the results of section V, where we have 
to resolve the scalar field as well.) As shown in |35| one 
can use such evolutions to obtain the frequencies of the 
excited modes of pulsation of the star by simply Fourier 
transforming those data. In general, however, the ex- 
citation of the different modes by the truncation error 
of the numerical schemes may not be sufficient to ac- 
curately determine the mode frequencies. Therefore, in 
order to compare unambiguously our mode frequencies 
with perturbation theory, we further perturb the density 
of the equilibrium models with an explicit eigenfunction 
p = p a +Ap c sin(7rr/_R), where R denotes the radius of the 
star (see Table I) and p is the density of the unperturbed 
star. The typical amplitude we use in the perturbation 
is A = 10~ 6 . 

Figures || and [| show the frequencies of the funda- 
mental mode and the first two overtones obtained by a 
Fourier transform of the radial velocity evolutions. The 
dashed vertical lines in these plots were obtained using 
a linear code |56| which evolves in time the perturba- 
tion equations. The agreement between the two codes is 
remarkable. The fundamental mode of the model with 
p c = 2.8 x 10~ 3 is already rather small, the star being 
close to the unstable branch. We note that the code is 
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FIG. 5. Fourier transform of the time evolution shown 
in Fig. |^. The frequencies of the fundamental radial mode 
(around / = 0.8 kHz) and the first two overtones are shown 
for a neutron star with central density p c = 2.8 x 10™ 3 . 

able to obtain a much higher number of overtones. Nev- 
ertheless, for the sake of clarity in the comparison and 
for our purpose of assessing the correct numerical im- 
plementation, it is sufficient to show only the first two 
harmonics. 



B. Scalar field dynamics in a regular spacetime 

In this section we present results aimed to validate the 
numerical implementation of the Einstcin-Klcin-Gordon 
solver. For this purpose we investigate the reflection of a 
scalar field at the origin of the coordinate system, turning 
off the hydrodynamics module of the code. The initial 
data for the scalar field packet are 

$ = 2xlO-V< r - 14 > 2 . (31) 

The location of this Gaussian pulse is chosen in such a 
way that, if superposed on the neutron star spacctimcs 
of the previous section, the scalar wave pulse would ini- 
tially lie outside the neutron star, the overlap being com- 
pletely negligible. Evolving this data, the initial pulse 
approaches the origin, is reflected, and radiates away, 
leaving behind Minkowski space. Such a sequence can be 
followed in Fig. |6[ for a simulation employing a grid of 
800 zones. We note the stability and smoothness of the 
solution, both at the origin and at J + . 

By evaluating global energy conservation, according to 
Eq. (p7[), after the pulse has reflected off the origin, we 



FIG. 6. Radial profiles of the scattering of a scalar field off 
the origin of the coordinate system at different times of its 
evolution. The numbers in the plot refer to the time coordi- 
nate u. For times u > 35 the pulse has completely radiated 
away leaving Minkowski spacetime behind. 

find that the energy is conserved (as expected) to second 
order accuracy. The best linear fit to the curve shown in 
Fig. [?] gives a slope of 1.99. 

As an aside we note that by simply changing the origin 
treatment in the code, it is possible to study the evolu- 
tion of a scalar field outside a spherical black hole. We 
performed such a simulation finding agreement with the 
results of @. 

C. Scalar field in a dynamic spacetime with a 
neutron star 

We consider now the full set of equations and prescribe 
initial data consisting of the ingoing scalar field pulse 
given by Eq. d3l|), together with a stable equilibrium, 
self-gravitating neutron star model with initial central 
density p c = 1.28 x 10" 3 , K = 100 and V = 2. This 
relativistic star model has a total mass of 1.4 Mq (see 
also p5fl ). We perform simulations of the scalar field scat- 
tering off the neutron star, focusing our study in this Sec- 
tion on the assessment of the global energy conservation 
properties of our complete numerical implementation. A 
comprehensive study of the dynamics of the scattering is 
deferred to Section |v|. 

Fig. |^ shows the Bondi mass of the neutron star - scalar 
field spacetime as a function of retarded time, combined 
with the total mass of the scalar field radiated away to 
null infinity. As one can clearly see from this figure, the 
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FIG. 7. Convergence test of the global energy conser- 
vation for the Einstein-Klein-Gordon system once the scalar 
wave has completely radiated away after reflecting off the ori- 
gin. The best linear fit to our results has a slope of 1.99, 
confirming the (global) second order accuracy of our numeri- 
cal implementation. 



spacetime is losing mass exactly at the rate which is ra- 
diated to null infinity by the scalar field. 

By computing Eq. ( |27| ) at a fixed retarded time of 
t = 0.5ms for different grid resolutions, we find that 
our code conserves globally the energy with a conver- 
gence rate which lies in between 1 and 2. The best linear 
fit is depicted in Fig. [|. The fact that the convergence 
rate drops now below second order is, however, to be ex- 
pected, since the approximate Riemann solver we are us- 
ing for the integration of the hydrodynamic equations is 
only (locally) first order accurate at discontinuities (i.e., 
the surface of the star) and at local extrema (i.e^, the cen- 
ter of the star) (see the related discussion in |||]). Nev- 
ertheless, for the highest resolution we have used, 2000 
radial grid points, the relative error in the energy conser- 
vation is of the order of 2 x 10 -6 for this very dynamical 
simulation. 




0.2 0.3 
retarded time [ms] 

FIG. 8. Bondi mass of the spacetime and total radiated 
energy of the scalar field as a function of retarded time for 
different resolutions. At the beginning, the scalar field con- 
tributes to the Bondi mass of the spacetime (1.4O5M0), before 
the Bondi mass drops in a small time interval, when the main 
part of scalar field mass (about 5 x 1O _3 M0) is emitted to 
future null infinity J + . As the sum of the two curves is con- 
stant, the energy is globally conserved (with a relative error 
of about 5 x 10~ 6 for the run with 1600 zones.) 



sian pulse of a scalar field according to Eq. thus 
fixing the pulse amplitude, width and location. For these 
initial data there is no overlap between the star and the 
scalar field pulse which then makes it possible to asso- 
ciate a specific initial mass with each one of the matter 
fields. 

In our setup we keep the central density as the only free 
parameter of the simulations, choosing a unique poly- 
tropic EoS and fixing the geometry and amplitude of the 
scalar field. This is clearly a severe restriction in the pa- 
rameter space of the scattering problem. Nevertheless, 
we choose this particular setup since we are interested 
in investigating the relativistic effects of the interaction, 
where the scalar field has a strong impact on the dynam- 



V. DYNAMICS OF SCALAR FIELD - NEUTRON 
STAR INTERACTIONS 

In this section we present our main results concerning 
the scattering of a scalar field pulse off relativistic neu- 
tron stars. As mentioned before, we use n = 1 relativistic 
polytropes as neutron star models. All models we con- 
struct lie on the stable branch of the total mass-central 
density diagram (see Fig. |l|) and are characterized by in- 
creasing central densities and compactness. Their basic 
properties are summarized in Table |]. 

In addition to the neutron stars we construct a Gaus- 



TABLE I. Equilibrium properties of the K = 100, n = I 
neutron star models in units in which c = G = Mq — 1. The 
entries are as follows: p c is the central density, M and R are 
the mass and radius of the star, respectively, and C = 2M/R 
is the compactness parameter. 



Pc (I0~ 3 ) M R C = 2M/R 

1.5 1.47 9.26 0.317 

2.2 1.60 8.45 0.379 

2.8 1.63 7.91 0.412 

2.9 1.64 7.84 0.418 
3.0 1.64 7.76 0.423 
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FIG. 9. Convergence test of the global energy conservation 
of a dynamic spacetime containing a self-gravitating neutron 
star and a scalar field. The rate of convergence is 1.62. See 
text for details. 



ics of the neutron star. A detailed analysis of the whole 
parameter space is beyond the scope of this work. 

When evolving in time the initial data, the scalar field 
travels inwards, enters the neutron star and it is finally 
reflected at the origin of the coordinate system. Con- 
trary to the Einstcin-Klcin-Gordon system wi thout the 
perfect fluid, which was discussed in Section IV B, the 



presence of the neutron star, and its associated potential 
well, may give rise to a phase of multiple interactions of 
the wave back and forth the origin and the maximum of 
the curvature potential. This, in turn, reflects itself in 
the existence of quasi-periodic signals (trapped modes), 
as discussed by [r2(| before its energy is radiated away. 
Furthermore, our neutron star models have been chosen 
conveniently close to the maximum of the stability curve 
(see Fig. [l]). Depending on the compactness of the neu- 
tron star onto which the wave pulse impacts, the stars 
are forced to either oscillate violently, or to collapse to a 
black hole on a dynamical timescale. Fig. [To] shows the 
spacetime diagram for the least compact neutron star 
model of our sample, with p c = 1.5 x 10~ 3 . For this 
model, the scalar field is able to force the star to con- 
tract and to expand, pulsating radially, as can clearly 
be identified in the varying location of the star's radius 
(the vertical solid line in Fig. |hJ). With our foliation, an 
outgoing scalar wave is covered in only one slice. 

Fig. O displays the time evolution of the central den- 
sity of the different neutron stars in our setup. The solid 
lines correspond to the neutron star-scalar field system. 
Correspondingly, the dashed horizontal lines indicate the 
evolutions of the equilibrium neutron star models without 
the presence of the scalar field. As already mentioned 
these evolutions are stable. They are characterized by 
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FIG. 10. Spacetime diagram of the reflection of a Gaus- 
sian scalar field pulse off a neutron star (K = 100, n = 1 and 
p c = 1.5 x 10~ 3 ). The dia gram focuses in the region close 
to the neutron star - scalar field interaction, but it was ob- 
tained from a global simulation of the spacetime. The dotted 
curves covering the whole diagram are outgoing light cones, 
which bend due to the spacetime curvature. The scalar field 
pulse, initially located at r — 14, travels inwards, enters the 
neutron star and is reflected at the origin of the coordinate 
system (the solid line corresponds to the maximum value of 
the scalar field). The interaction with the scalar field triggers 
the oscillation of the neutron star, which can be seen from the 
vertical solid line of varying location in the diagram, which 
indicates the radius of the star. 



the appearance of small-amplitude oscillations associated 
with the radial modes of pulsation of the star (which are 
too small to be seen in the figure). On the other hand, 
all neutron star-scalar field models with initial central 
density below 2.8 x 10 -3 also oscillate around the stable 
equilibrium model. The oscillation frequencies of the two 
least compact models correspond to the frequencies cal- 
culated in the linear regime, even though the amplitude 
of the oscillations is now much larger due to the scalar 
field kick. This is no longer the case for the model with 
a central density of p c — 2.8 x 10~ 3 . For this model, 
which is close to the threshold of black hole formation, 
the amplitude of the oscillations is big enough to show 
nonlinear effects, the oscillation frequency being much 
smaller than the value obtained in the previous section. 
For the models with central densities of p c = 2.9 x 1CP 3 
and p c = 3.0 x 10 -3 the interaction with the scalar field 
is able to trigger their gravitational collapse to a black 
hole on a dynamical timescale. Unfortunately, due to nu- 
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FIG. 11. Central density of the neutron stars interacting 
with the scalar field as a function of retarded time. The 
three less compact models, with p c < 2.8 x 10~ 3 , oscillate 
strongly around their equilibrium value after interacting with 
the scalar field. The other two more compact models collapse 
to a black hole instead on a dynamical timescale. The dashed 
lines are taken from our evolutions of the equilibrium model 
without the presence of the scalar field. 

merical inaccuracies arising at the end of the simulation 
we are not able to follow the collapse process once the 
event horizon is about to form (similar problems were 
reported in |27j for the collapse of supermassive stars). 
However, convergence studies show clear evidence that 
these models collapse to black holes. Further evidence is 
given by the evolution of the neutron star radii, as shown 
in Fig. [12|. 

To demonstrate the dynamic range of these evolutions, 
we focus on the model with central density p c — 2.9 x 
1CP 3 . Initially, the redshift factor e 2H between the center 
of the star and observers located at r ~ * oo is 2.1. By 
the end of the simulation it has increased to a value of 
59.5. This high redshift factor explains the appearance 
of a kink in the central density towards the end of our 
numerical evolution (see Fig. |l l| ). 

We also note that global energy conservation is very 
well fulfilled for these extreme hydrodynamic simula- 
tions. The relative deviation from energy conservation 
according to Eq. ( p7| ) when the evolution stops is of the 
order of 10 -4 . 

By analyzing the energy transfer from the scalar field 
to the neutron star during the interaction we find that it 
increases the more compact the stellar models are. This 
behavior is shown in Table [Tl| We remark that the initial 
mass of the scalar field is not strictly identical for all cases 



FIG. 12. Time evolution of the radius of the different neu- 
tron stars interacting with the scalar field. The radius of two 
most compact models decreases dramatically, indicating that 
they undergo gravitational collapse to a black hole. 

considered, due to the different underlying geometry on 
which the Gaussian pulse is constructed. We evaluate 
the total radiated mass in the scalar field at a retarded 
time of r = 0.6ms. The mass radiated away to infinity 
after this time is negligible. 

Next we analyze the behavior of the scalar field in these 
scattering simulations. In Fig. Il3| w e plot the (retarded) 
time evolution of the news, Eq. (|26|) , for the whole sam- 
ple of our neutron star models. The scalar field signal 
measured at null infinity can be neatly divided into three 
phases. The first phase, before the main pulse reflects 
off the origin (not shown in the figure), is dominated by 
an initial backscattering, the amplitude of the signal be- 
ing small. The second phase, whose duration depends on 
the compactness of the neutron star |26| , is characterized 

TABLE II. Energy transfer from the scalar field to the 
neutron star during the scattering process. The entries are as 
follows: p c is the central density of the neutron star, M* is 
the initial mass of the scalar field, £' rac j is the total radiated 
mass, and £j rans is the percentage of the energy transfered 
in the interaction. We use units in which G = c = Mq = 1. 



Pc (10- 3 ) 


M * (10- 3 ) 


^rad (W- 3 ) 


^trans 


1.5 


4.90 


4.86 


0.8 


2.2 


4.80 


4.72 


1.7 


2.8 


4.76 


4.65 


2.3 


2.9 


4.75 


4.63 


2.5 


3.0 


4.75 


4.62 


2.7 
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by the reflection of the main scalar field pulse back and 
forth the origin and the maximum of the neutron star 
curvature potential, which, in turn, induces the appear- 
ance of quasi-normal oscillations on the scalar field. Most 
of the energy is radiated away in this period. Once the 
pulse has lost sufficient energy it enters a third phase, in 
which the late time behavior of the signal is dominated 
by a power-law tail N oc t~ a , with a = 3, due to the 
reflection of the scalar field at the exterior Schwarzschild 
geometry pT| , |i7| , |l8| , p6| . Since the compactness of our 
models is well below the Buchdahl limit, C — 8/9, the 
quasi-normal mode ringing phase does not last for an 
extended period of time. Therefore, after a few reflec- 
tions trapped inside the curvature potential, the signal 
enters rapidly the power-law tail phase. From Fig. ^ 
one can see that the more compact the neutron star, the 
larger the quasi-normal ringdown phase. We also point 
out that by going to more compact models, increasing 
the central density of the neutron star beyond the max- 
imum of the stability curve (i.e., going into the unstable 
branch) and freezing the hydrodynamics and metric evo- 
lution to avoid gravitational collapse, we are able to find 
a much longer ringdown phase. Our results, obtained for 
fully self-gravitating, polytropic neutron star models, are 
in perfect agreement with previous findings by Pavlidou 
et al [Bq] , who used a more idealized setup consisting of 
constant density, static neutron stars. 

Since the study of the late time power-law tails re- 
quires sufficient resolution, especially for large radii, we 
have used a different radial coordinate for these simu- 
lations, r = 3Cte/(l — x 4 ). This allowed us to resolve 
the power-law behavior, avoiding the evolution from be- 
ing dominated by numerical noise. By performing a 
linear regression study of the tails in the time interval 
log(r \ms\) € [0.3; 0.7], we obtain the results summarized 
in Table III. We find the correct power-law behavior 



of the scalar field in our fully dynamical evolutions, as 
predicted by both, linear analysis and by nonlinear nu- 
merical evolutions of scalar fields in the exterior black 
hole geometry |l9|| - Note that we measure the tails 
on the news, whereas the results of the above references 
read off the quantity g at future null infinity . Both 
quantities are related by 
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FIG. 13. Time evolution of the news function during the 
scattering problem. The different lines correspond to the dif- 
ferent models in our sample of Table |, and are labeled in the 
plot with respect to the compactness parameter. The dura- 
tion of the more dynamic quasi-normal ringing phase strongly 
depends on the compactness of the neutron star model, in- 
creasing as the compactness increases. The late time behavior 
of the signal decays as an inverse power-law. 



studied the fate of the neutron stars when they are hit 
by a sufficiently strong scalar field packet, as well as the 
dynamics and energetics of the process. 

We have found that by choosing a strong (finite am- 
plitude) scalar field pulse with a mass of the order of 
10~ 3 Mq, the neutron star is either forced to oscillate in 
its radial modes of pulsations or to collapse to a black 
hole on a dynamical timescale. The fate of the star de- 
pends on its central density and, since we fix the poly- 
tropic equation of state, on its compactness. The energy 
transfered to the neutron star increases with the com- 
pactness of the model. The radiative signals we have 
found in our fully nonlinear simulations consist of sev- 



VI. SUMMARY 

We have analyzed numerically the interaction of neu- 
tron stars and scalar fields by means of nonlinear evolu- 
tions of the Einstein-Klein- Gordon perfect fluid system in 
spherical symmetry. We have built a sequence of stable, 
self- gravitating, K = 100, n = 1 relativistic polytropes, 
increasing the central density from p c = 1.5 x 10~ 3 to 
3.0 x 10~ 3 (G = c = Mq = 1). Using a compacti- 
fied spacetime foliation with outgoing null cones we have 



TABLE III. Late time 
power-law behavior of the news 
N cc t~ a for the (stable) neu- 
tron star - scalar field scattering 
problem. The results agree with 
the predicted exact value a — 3. 
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eral quasi-normal oscillations and a late time power-law 
tail, in agreement with the results predicted by (linear) 
perturbation analysis of wave propagation in an exterior 
Schwarzschild geometry |21[. 
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